Using R in hydrology (EGU2017 short course)

These slides and all other course materials can be found at
github.com/brry/rhydro #slides
To get all the material including the datasets and presentation source code, we recommend to download the whole github course repository.

This is an R Markdown Notebook.
For discussions, please visit the Hydrology in R Facebook group.
Before running the code blocks below, we suggest to get package installation instructions by running:

source("https://raw.githubusercontent.com/brry/rhydro/master/checkpc.R")

Aim and contents of this workshop

We want to:

• Show off how awesome R is for hydrology (it’s R-some!^^)
  
• Convince you to start or continue using R
  
• Provide all the code for you as a starting point

We can not:

• Teach you actual R coding (90 mins is too short for a tutorial)

We have prepared:

• Good coding practice, report generation (Rstudio, rmarkdown, R notebook)
• Using R as GIS (reading a rainfall shapefile + Kriging, sf + leaflet + mapview + OSMscale)
• River discharge time-series visualisation and extreme value statistics (animation + extremeStat)
• Hydrological modelling with airGR
• Exploratory Data Analysis including flow duration curve and trend analysis on time-series

 

Before we get started, please let us know your current R knowledge level by filling out the short survey at
bit.ly/knowR

top

Report

Good coding practice, report generation (Rstudio, rmarkdown, R notebook)
Daniel Klotz

top

GIS

Using R as GIS (reading a rainfall shapefile + Kriging, sf + leaflet + mapview + OSMscale)
Berry Boessenkool

Shapefiles

Reading shapefiles with maptools::readShapeSpatial and rgdal::readOGR is obsolete.
Instead, use sf::st_read. sf is on CRAN since oct 2016.
Main reaction when using sf: “Wow, that is fast!”
Download the shapefile or better: download the whole github course repository

rain <- sf::st_read("data/PrecBrandenburg/niederschlag.shp")

Reading layer `niederschlag' from data source `/home/berry/Dropbox/R/rhydro/presentations/data/PrecBrandenburg/niederschlag.shp' using driver `ESRI Shapefile'
converted into: POLYGON
Simple feature collection with 277 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 3250175 ymin: 5690642 xmax: 3483631 ymax: 5932731
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=3500000 +y_0=0 +ellps=GRS80 +units=m +no_defs

centroids <- sf::st_centroid(rain)
centroids <- sf::st_coordinates(centroids)

top

Plotting, maps

Static plot:

plot(rain[,1])

Static map:

prj <- sf::st_crs(rain)$proj4string
centroids <- as.data.frame(centroids)
cent_ll <- OSMscale::projectPoints(Y,X, data=centroids, to=OSMscale::pll(), from=prj)
map_static <- OSMscale::pointsMap(y,x, cent_ll, fx=0.08, type="maptoolkit-topo", zoom=6)

Downloading map with extend 11.029332, 14.981464, 51.10418, 53.765423 ...
Done. Now plotting...

Interactive map:

library(leaflet)
cent_ll$info <- paste0(sample(letters,nrow(cent_ll),TRUE), ", ", round(cent_ll$x,2), 
                       ", ", round(cent_ll$y,2))
leaflet(cent_ll) %>% addTiles() %>% addCircleMarkers(lng=~x, lat=~y, popup=~info)

Interactive map of shapefile:

# devtools::install_github("environmentalinformatics-marburg/mapview", ref = "develop")
library(berryFunctions) # classify, seqPal
col <- seqPal(n=100, colors=c("red","yellow","blue"))[classify(rain$P1)$index]
mapview::mapview(rain, col.regions=col)

top

Kriging

Plot original points colored by third dimension:

pcol <- colorRampPalette(c("red","yellow","blue"))(50)
x <- centroids$X # use cent_ll$x for projected data
y <- centroids$Y
berryFunctions::colPoints(x, y, rain$P1, add=FALSE, col=pcol)

Calculate the variogram and fit a semivariance curve

library(geoR)

--------------------------------------------------------------
 Analysis of Geostatistical Data
 For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
 geoR version 1.7-5.2 (built on 2016-05-02) is now loaded
--------------------------------------------------------------

geoprec <- as.geodata(cbind(x,y,rain$P1))
vario <- variog(geoprec, max.dist=130000) # other maxdist for lat-lon data

variog: computing omnidirectional variogram

fit <- variofit(vario)

variofit: covariance model used is matern 
variofit: weights used: npairs 
variofit: minimisation function used: optim 

initial values not provided - running the default search

variofit: searching for best initial value ... selected values:
              sigmasq   phi        tausq kappa
initial.value "1325.81" "19999.05" "0"   "0.5"
status        "est"     "est"      "est" "fix"
loss value: 104819060.429841 

plot(vario)
lines(fit)

Determine a useful resolution (keep in mind that computing time rises exponentially with grid size)

# distance to closest other point:
d <- sapply(1:length(x), function(i)
            min(berryFunctions::distance(x[i], y[i], x[-i], y[-i])) )
# for lat-long data use (2017-Apr only available in github version of OSMscale)
# d <- OSMscale::maxEarthDist(y,x, data=cent_ll, fun=min)
hist(d/1000, breaks=20, main="distance to closest gauge [km]")

mean(d/1000) # 8 km

[1] 8.165713

Perform kriging on a grid with that resolution

res <- 1000 # 1 km, since stations are 8 km apart on average
grid <- expand.grid(seq(min(x),max(x),res),
                    seq(min(y),max(y),res))
krico <- krige.control(type.krige="OK", obj.model=fit)
#krobj <- krige.conv(geoprec, loc=grid, krige=krico)
#save(krobj, file="data/krobj.Rdata")
load("data/krobj.Rdata") # line above is too slow for recreation each time

Plot the interpolated values with or an equivalent (see Rclick 4.15) and add contour lines.

par(mar=c(0,3,0,3))
geoR:::image.kriging(krobj, col=pcol)
colPoints(x, y, rain$P1, col=pcol, legargs=list(horiz=F, title="Prec",y1=0.1,x1=0.9))

points(x,y)
plot(rain, col=NA, add=TRUE)

top

Discharge

River discharge time-series visualisation and extreme value statistics (animation + extremeStat)
Berry Boessenkool

top

Hydmod

Hydrological modelling with airGR
Katie Smith

top

EDA

Exploratory Data Analysis including flow duration curve and trend analysis on time-series
Shaun Harrigan

top

Discussion

Please give us feedback at bit.ly/feedbackR

For discussions, please use the Hydrology in R Facebook group.